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Abstract. - The effect of coordination on transport is investigated theoretically using random 
networks of springs as model systems. An effective medium approximation is made to compute 
the density of states of the vibrational modes, their energy diffusivity (a spectral measure of 
transport) and their spatial correlations as the network coordination z is varied. Critical behaviors 
are obtained as z — > z c where these networks lose rigidity. A sharp cross-over from a regime where 
modes are plane-wave-like toward a regime of extended but strongly-scattered modes occurs at 
some frequency uj* ~ z — z c , which does not correspond to the Ioffe-Regel criterion. Above uj* 
both the density of states and the diffusivity are nearly constant. These results agree remarkably 
with recent numerical observations of repulsive particles near the jamming threshold [T|. The 
analysis further predicts that the length scale characterizing the correlation of displacements of 
the scattered modes decays as l/y/u with frequency, whereas for uj « uj* Rayleigh scattering 
is found with a scattering length l s ~ (z — z c ) 3 /uj 4 . It is argued that this description applies to 
silica glass where it compares well with thermal conductivity data, and to transverse ultrasound 
propagation in granular matter. 



The temperature dependence of the thermal conductiv- 
ity in glasses above a few K° suggests a regime where the 
scattering of phonons is strong and nearly independent of 
frequency, as noticed by Kittel [5] . At these temperatures 
specific heat measurements indicate the presence of excess 
modes in comparison with the Debye model, the Boson 
Peak [3], whose amplitude correlates with other funda- 
mental properties of glasses, such as the fragility of the 
super-cooled liquid [J. Understanding the nature of the 
strongly scattered modes leading to an excess density of 
states may therefore unravel some of the most mysterious 
properties of these materials. Several mean-field models 
have been proposed to explain the presence of such modes 
[5HI] and their effects on transport. Nevertheless, there 
is yet no simple generic rules to predict from the micro- 
scopic structure of the glass the amplitude of the peak 
and the frequency at which transport becomes diffusive. 
Existing mean-field approaches do not explain why, for ex- 
ample, silica glass presents one of the largest Boson peak, 
whereas it is almost inexistent for amorphous silicon, al- 
though both structure are tetravalcnt. The spatial nature 



of the modes also remains unresolved. Recently numer- 
ics have shown that the vibrational properties of repul- 
sive, short-range particles display scaling behavior when 
decompressed toward the jamming transition [5]. One ob- 
serves a frequency scale uj* which vanishes with decom- 
pression and marks a sharp cross-over toward a regime of 
strong scattering at larger frequency, where the modes dif- 
fusivity, which can be computed numerically via a Kubo 
formalism [TJ[H], is nearly constant and small, in consis- 
tence with Kittel's conjecture, and where the density of 
states displays a plateau. Close to the jamming thresh- 
old the cross-over toward strong scattering can occur at a 
wave-length infinitely larger than the particle size. Thanks 
to this separation of length scales this system is very suit- 
able to study the spatial structure of the excess mode, and 
its dependence on the microscopic structure. A variational 
argument can explain the plateau in the vibrational spec- 
trum, and shows that both the mean coordination and 
contact strain, averaged on a scale I* ~ l/(z — z c ) which 
diverges near threshold, control locally the onset of the 
excess modes [TUHH]- referred to as "anomalous" to un- 
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derline the difference of these excitations with plane waves. 
This length does not appear in the static structure, but 
characterizes the heterogeneity of the response to a point 
force [12] • Nevertheless, this variational argument does 
not yield the spatial correlations of the displacements of 
the modes, nor their transport ability, or diffusivity. Inter- 
estingly, some aspects of the diffusivity can be estimated 
by making several assumptions on the spatial nature of 
the modes or their wave- vector decomposition [9]. How- 
ever, these assumptions are not trivial and lack theoreti- 
cal justifications. In this letter the spatial correlations of 
the modes and their diffusivity are computed by applying 
an effective medium approximation, which has been used 
in similar contexts [TUGS] and in electrical resistor net- 
works [16j , to random networks of harmonic springs close 
to the minimal coordination where rigidity is lost and un- 
jamming occurs. Results compare well with properties of 
sphere packing, and also with measurements of thermal 
conductivity in silica at high temperature. 

I consider a system of point particles of unit mass 
m = 1, interacting via a lattice of springs with coor- 
dination zq strictly larger than twice the spatial dimen- 
sion d, where all bonds are assumed to be equivalent, 
such as in the triangular lattice. Springs are harmonic, 
unstretched, of unit length a = 1, and have a stiff- 
ness ko, whose square root defines our unit frequency 
(thus fco = 1 in our units). The energy expansion reads 

SE = Yl/ij) |[0^> - <^j) ' ™ij\ 2 + o(5R 2 ), where the sum 
is on all bonds ij, whose direction is along the unit vec- 
tor Hij . Springs are randomly removed with a probability 
1 — p, such that the coordination of the remaining net- 
work is z = pzQ. The static clastic properties and the 
presence of vibrational modes of zero frequency of such a 
disordered lattice have been studied previously [2]. Here 
we extend the analysis to compute vibrational properties 
in the mechanically stable phase. 

The random network is modeled as an effective medium 
of spring stiffness fcj\/(w). To determine fcA/(cj), one con- 
siders the "defectuous" network where all bonds have a 
stiffness /cm(w), except for one bond (ij) = (12) at the 
origin, whose stiffness fci2 is distributed according to the 
distribution P(x) = pS(x = 1) + (1 — p)5(x = 0) of the dis- 
ordered system. Next one considers a homogeneous oscil- 
latory stress applied on this system. If the bond (12) were 
identical to the other bonds, the applied stress would lead 
to an affine oscillatory response 8Rf(t) = 5Rf(oj) exp(zwt), 
i = 1...N . Such a response is not a solution in general, as 
if it is imposed on this network, an unbalanced oscillatory 
force dipole of amplitude fi 2 {ui) = (fci2 — k]\,{(uj))e(uj) ap- 
pears on the bond 12, where e(oj) — (5Ri(tju)—5R,2(w))-n\2 
is the contact strain. The solution is obtained by adding 
the response of this dipole to the affine displacement. This 
causes an extra strain e*(w) to appear in the bond (12). 
/cm(w) is determined self-consistently by imposing that the 
effective medium respond with the correct average strain 
to an imposed stress, implying that (e*(to)) = 0, where 



(X) = J X(ki2)P(ki2)dki2 . In general this requirement 
leads to a complex solution for /cm(w) , reflecting the bro- 
ken symmetry of translation in the disordered network and 
the presence of scattering. 

The extra strain e*(w) can be expressed in terms of 
the Green function of the effective homogeneous system 
[TUHnjITT], as is now sketched. Newton's equation for an 
oscillatory external force field SFi(t) = 5Fi(uj) exp(iwi) 
leads to J2j ^y( CJ ) ' ( ^^i( a; ) = SFi(uj), where: 

Aij(iL>) = 8 lj (~u 2 + ^kMnii®riu)~5(ij)k M riij®ni : j (1) 
(i) 

In Eq.(UJ ® is the tensor product, the sum is made on all 
neighbors / of i, and 8uj\ = 1 if i and j are neighbors, and 
zero otherwise. In these notations Aij is a tensor of rank 
d. We denote by A(uj) the dN x dN linear application 
relating displacements to forces field. Its inverse is the 
Green function G = A -1 . Gij is the tensor of rank d 
characterizing the displacement of particle i in response 
to a force applied in j. In the "defectuous" network the 
additional displacement field 6R*, i — 1...N, generated by 
the force dipole in the bond 12 satisfies Vi: 

Y^AijSRj - (hi - 5 l2 ){k 12 - fc M )ni2 ® n 12 ■ (SR\ - 6R 2 ) 
j 

= -ni2{k M {w) - ki 2 )e(io)(5ii - 5 i2 ) 

After some algebra one gets: 

e*{u) = {5R\{u)~5R* 2 {Lu))-n 12 

2e(tL>)(k M - k 12 )n 12 ■ [Gujuj) - Gi 2 {w)) ■ n 12 
1 + 2(ki 2 - k M )ni 2 ■ (Gii(w) - G 12 (u>)) ■ n 12 

This expression can be simplified using that AijGjiis 
the unit tensor of rank d, Eq.(UJ), and the assumption that 
all bonds are equivalent [16]. One obtains if i and j are 
neighbors that «y • (Gh(lo) — Gi 3 (w)) • n,j = d(l + riij ■ 
Gn(bj) ■ Hij) I '(kMZo). Applied in Eq.(J5]), and requiring 
(eo(w)) = leads to: 



k M + (ki2 - fc M )|f(l + w 2 ni2 • G n (a;) ■ n 12 ) 

Note that in this expression, n\2 can be replaced by any 
contact direction Hij and Gn(w) by Gu(iu). 

Taking the static limit uj — > in Eq.((2]) and using pzo = 
z, one gets k^[ = (z — 2d)/(z — 2dp). Thus as found 
by Thorpe and collaborators [33], the effective stiffness of 
the network, and therefore all the elastic constants, vanish 
linearly as z — > 2d, as found numerically in generic elastic 
networks [15] . z = 2d corresponds to the Maxwell bound 
where rigidity is lost [18], such network are called isostatic. 
Note that for sphere packing, the bulk modulus does not 
vanish near the isostatic point, as explained in [3Tj. For 
such purely repulsive systems I expect the forthcoming 
result on the Rayleigh scattering to hold for transverse 
modes only. 
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The Green function can be computed numerically for 
any lattice, and used to solve the self-consistent Eq.©. In 
what follows we focus on the low-frequency scaling prop- 
erties near this critical point. One expects such properties 
to be universal and not to depend on the details of the lat- 
tice. We shall therefore consider a simplified model, where 
the lattice elasticity is characterized by a single clastic 
constant B, which is taken equal to km- For the isotropic 
lattices, such as the hexagonal lattice, two elastic constant 
are required to describe elasticity. The analysis below is 
readily extended to this case, and leads to the same scaling 
behavior. We consider for the Green function: 



Gii(w) 



d d k 
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fcA/(w)fc 2 — UJ- 



;k a k a (3) 



where the Debye cut-off k max is a constant in our units, 
and is set to one in what follows, r,; is the vector going 
from particle 1 to i, and the sum is made on the different 
polarization k a . All spatial directions being equivalent, 
one obtains: 



trace((?ii (u>))/d = H12 ■ Gu(uj) ■ Hi 
d d k 1 



o (2u) d k M (uj)k 2 - uj 2 



(4) 
(5) 



Eas. ([2l5p are solved numerically in the three-dimensional 
case, using z = 12. The results for the real and imag- 
inary part of &m = ki(uj) — iY,(ui) are shown in Fig.([T]) 
for three different coordination. Qualitatively these re- 
sults can be obtained by approximating the real part of 
the integral of Eq.© by its large k contribution (obtained 
by setting uj = 0), and its imaginary part by the singular 
value of the integral, plus its large k contribution. This is 
approximately valid for k\ > S. Expanding Eq. (|2J) at low 
frequency, one finds the following schematic set of equa- 
tions (where numerical pre-factors are omitted): 



ki = 8z 



3/2 



UJ 2 kl 

£ 2 



k\ 



(6) 
(7) 



where 6z = z — z c characterizes the distance from thresh- 
old. These schematic equations lead to a sharp cross-over 
at uj* = Sz/2. This scaling is checked in inset of Fig.©. 
For smaller frequencies, k\ sa Sz, as expected from the 
static analysis, and E m uj 3 /5z 3 / 2 « ki. At u*, 
drops by some numerical factor, and S increases by a fac- 
tor of order l/y/6z. For larger frequencies, £ >> ki, and 
the (middle) term of Eq. © corresponding to the singular 
value of the integral is not present. This leads to S ~ to. 
The density of states D(uj) of the effective medium is 
readily extracted using D(uj) = (2w/7r)Im(tracc(Gn(a;)), 
and is shown in Fig.©. The cross over at uj* cor- 
responds to a sharp transition from a Debye behavior 
D(uj) ~ <jz~ 3 l 2 uj 2 at low frequency toward a plateau at 
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Fig. 1: Logarithmic plot of the real and imaginary part of 
the effective stiffness ku = ki — iE vs frequency for several 
coordination numbers Sz = z — z c . 



higher frequency, and therefore correspond to the boson 
peak. The same cross-over is observed numerically for 
elastic particles near jamming 

We next consider the spatial Green function at fre- 
quency uj. For large distances fj >> 1 one finds: 



trace(Gii(w)) 



Anvri k 



(8) 



I s (oj) and v(u>) can be viewed as the scattering length and 
the velocity of the corresponding vibrations, and their ex- 
pression is: 



v(oj) 



oj\lm(k]^) 

\k M \ 
Rc(fclf) 



(9) 
(10) 



Transport at frequency uj is characterized by the diffu- 
sivity Di(u), whose approximation Di(uj) = v(uj)I s (uj)/3 
is shown in FigQ. This approximation is expected to 
be accurate for weak-scattering (uj << uj*) and crude 
but qualitatively correct for the strongly-scattered modes 
(uj » uj*) as long as they are extended, since it enforces 
the scaling requirement that the diffusivity is the product 
of a characteristic speed times a characteristic length. A 
more accurate computation of the diffusivity would require 
to go beyond the single bond effective medium investigated 
here. 

For uj « uj* . Rayleigh scattering occurs: Di(uj) ^ 
6z 7 / 2 /uj 4: . At uj uj* , there is a sudden decay of the 
scattering length from 1/Sz toward 1/s/Sz which is of the 
order of the wavelength A = v/uj. This sudden decay of the 
scattering length at uj* implies that this cross-over cannot 
be identified from the quasi-plane wave regime (uj << uj*) 
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log(ffl) 

Fig. 2: Density of states vs frequency uj for five coordination 
numbers. Inset: Log-log plot of ui* , defined as D(ui*) = 0.2, vs 
Sz. The line has a slope 1. 

using the Ioffc-Regel criterion l a = A. For uj >> uj* there 
is a flat plateau in the diffusivity, again in very good agree- 
ment with numerical studies of transport of jammed elas- 
tic particles [2(2]. This plateau stems from the compen- 
sation of two effects: the decrease of the scattering length 
l s (uj) ~ cj -1 / 2 , and the inverse increase of the velocity 
v(uj) ~ uj 1 / 2 . Note that at the high frequency end of the 
spectrum a localization edge, rather insensitive to coordi- 
nation, is observed numerically [T] and is not captured by 
this analysis. 

The spatial properties of the strongly scattered normal 
modes \e u ) are deduced from the relation (e u (l)(8)e tJ (i)) cx 
Iui(G u {uj))uj/D(uj). Thus Eq.© implies that at uj*, the 
correlation length of the displacement is of order <5z -1 / 2 , 
and decreases as ui~ x l 2 for higher frequencies. Scaling ar- 
guments for why anomalous modes near uj* are character- 
ized by a length of order Sz^ 1 / 2 -which is the value of the 
wavelength just below uj* [T2j- have been given by Vitelli 
ct al. 0. 

Experimentally the dynamical structure factor S{k, uj) 
is accessible. In Fig.(j4|) the quantity Im(A: 2 / (k^jk 2 — w 2 )), 
proportional to ujS(k, uj) for a harmonic dynamics, is plot- 
ted. At threshold, no peaks are found in the structure 
factor, which simply consists of a gap at low frequency, 
which ends approximately when uj 4 ~ E 2 fc 4 leading to 
k ~ uj 1 / 2 , before reaching a plateau. At larger coordina- 
tion, well defined-peaks for uj < uj* appear in the spec- 
trum, corresponding to weakly-scattered plane waves. At 
larger frequencies the peak erodes, and the structure fac- 
tor eventually resembles what is found at the critical point. 
The intensity gap at small frequency and the presence of 
a large background of large spatial frequencies have been 
observed above the Boson peak peak in various systems 
such as silica [20] and sphere packing [15] . 

I define a dispersion relation by considering the maxi- 
mum k max (uj) of S(k, uj) at fixed uj. As shown in the lower 
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Fig. 3: Log- log plot of the diffusivity Dt(ui) vs ui for several 
coordination numbers, as indicated in the legend. 

panel of Fig. ((3]), this dispersion relations is linear at low 
frequency, and presents an inflection point near the Boson 
peak, before increasing as uj ~ k 2 . This inflection point 
leads to a drop in the velocity uj / k max at the Boson peak 
frequency as shown in the inset of Fig.®, a rather subtle 
effect recently observed in Lennard- Jones simulations [H| . 

Applications to amorphous solids: Numerical observa- 
tions show that when friction is added, the scaling rela- 
tions valid for the vibrational spectrum of purely repulsive 
particles still hold [22]. Thus the results on transport de- 
rived here are expected to hold in granular matter as well, 
at least in the limit where dissipative effects [23] do not 
dominate sound propagation. In small systems of grains, 
low-frequency ultrasounds propagate ballistically, whereas 
propagation is diffusive at high-frequency [2_3] . According 
to the present analysis the cross-over wavelength between 
these regimes should decrease with coordination. The co- 
ordination depends on the system preparation and the 
friction coefficient, and can be measured independently 
by comparing transverse and longitudinal speed of sound 
[25H27]. 

Covalent glasses are also well suited to study the effect 
of coordination on transport. Practically important ex- 
amples are tetravalent covalent networks where the joints 
linking tetrahedra have a small bending force constant 
such as in silica glass (where the Si-O-Si angle is eas- 
ily deformable) or germanium oxide. If these joints are 
are assumed to be freely rotating and only the other co- 
valent interactions are taken into account (the so-called 
RUM model [25] ) , such networks turns out to be isostatic 
and display a flat density of states up to zero frequency 
|29j . in consistence with the present analysis and previ- 
ous arguments [301I31] . When the small but finite bending 
force constant among tetrahedra is taken into account, the 
plateau of anomalous modes in the spectrum is shifted and 
appears above IThz [3D1I2I], corresponding to the Boson 
peak frequency. Most of the spectrum thus consists of 
anomalous modes which must display a constant diffusiv- 
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Fig. 4: Dynamical structure factor S(k, uj) times u vs wave 
vector k for different frequencies indicated in legend. Upper 
panel: Sz = 0.29. Middle panel: Sz = 4 * 1CT 4 . No peaks 
are apparent. Lower panel: Dispersion relation fcmax(w) as 
defined in the text for Sz — 0.29. Inset: phase velocity oj/k ma x 
vs frequency w for Sz = 0.29. The velocity displays a minimum 
close to the Boson peak frequency to bp = 0.09. 

ity according to the present analysis. Assuming a constant 
density of states and diffusivity up to some maximal fre- 
quency uj m ax enables to compute the thermal conductivity 
K ( T ) = f™ D(u)C(u)Di(u)du oc_ J^ mB " C(w)doj where 
C(lu) oc (a;/T) 2 exp(/iw/fcT)/(exp(/ia;/fcT)-l) 2 istheheat 
capacity of the modes at frequency u>. Fig.© shows that 
such an estimate captures accurately the thermal conduc- 
tivity for several decades above its plateau. In this high- 
temperature estimate, I have neglected the contribution 
of plane- wave like modes at low frequency. This is a good 
approximation since the contribution of these modes can 
be read at lower temperature where only those modes af- 
fect transport and is small. The plane waves contribution 
is of order of the thermal conductivity plateau where the 
cross-over occurs, which is significantly smaller than the 
high temperature conductivity. This analysis thus sup- 
ports that anomalous modes are responsible for transport 
at room temperature. Note that in tetravalent structures 
with strong bending force constant at the joints, such as 
in amorphous silicon where the same covalent bonds char- 
acterize the joints and the tetrahedra themselves, the ef- 
fective coordination is high, and a Boson peak is hardly 
detectable, as expected. The present results are expected 



Fig. 5: Log-log plot of the thermal conductivity in silica glass 
against temperature in arbitrary units. Observations taken 
from Freeman and Anderson [33] appear in crosses. Full line: 
prediction for a flat spectrum of anomalous modes with con- 
stant diffusivity. 

to apply in simple molecular liquids where long-range in- 
teractions are present as well |3"Tll32) . 

Comparison with the variational argument of [10]: In 
|10] a variational argument guessing the nature of the 
modes leading to a plateau in the density of states was 
able to explain the abrupt change in this quantity at some 
frequency uj* ~ 5z. The corresponding anomalous modes 
are not built by softly distorting translational modes like 
plane waves are. The view that the Boson peak is a cross- 
over between two distinct types of excitations is reflected 
by the sharp cross-overs in the density of states and the 
transport properties at ui* found here, supporting that 
new objects suddenly appear in the spectrum at this fre- 
quency. Although anomalous modes also can be charac- 
terized by a length scale and a velocity, their spatial struc- 
ture qualitatively differs from plane waves as appears in 
their dynamical structure factor, flat with a gap at small 
frequencies rather than peak-like as for plane waves. 

If the variational procedure of [TO] and the effective 
medium theory developed here lead to consistent results, 
these approaches are also complementary. The variational 
argument suggested that anomalous modes would lead to 
poor transport [9ll31j , but did not enable to compute their 
dynamical structure factor explicitly nor the correlation 
length l s ~ l/y/Sz of their displacement field, as was done 
here. One the other hand, the present approach does not 
yield the length scale /* ~ 1/Sz on which fluctuations in 
the microscopic structure affect the anomalous modes [TU] 
and the response to a local perturbation [12]. This may 
be a generic limitation of mean field approaches. 

Comparison with previous mean field approaches: Pre- 
vious mean- field analysis of the Boson peak, see e.g. [SHI], 
have also connected this phenomenon to the presence of 
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an elastic instability. Nevertheless the nature of the insta- 
bility in those cases is singularly different from the present 
study, for example they find that the density of states does 
vanish at zero frequency at the instability threshold. This 
discrepancy underlines that two parameters, and two fre- 
quency scales, must be included to describe excess modes 
even in simple amorphous solids. In particular, the present 
analysis has not included the effect of compression on the 
vibrational spectrum, or equivalently the presence of a fi- 
nite strain in the contacts. This is a good approximation 
for covalent glasses and granular matter at limited pres- 
sures, but compression plays an important role e.g. for col- 
loids near the glass transition where it even determines the 
microscopic structure of the glass [21] ■ For any fixed coor- 
dination z > 2d, increasing compression shifts anomalous 
modes toward lower frequencies and eventually leads to an 
clastic instability as shown in for which the density 
of states does vanish at zero frequency. Near this insta- 
bility, the spectrum displays two decoupled characteristic 
frequency: a Boson peak frequency where the first anoma- 
lous modes appear and a higher frequency w* above which 
the density of states displays a plateau [TT] , and where the 
present analysis is expected to hold. These two character- 
istic frequencies are identical only when pressure is negli- 
gible, in general they can be different, reflecting that two 
parameters, coordination and pressure, strongly affect the 
vibrational spectrum. Previous mean-field analysis [5 7 
are presumably well suited to describe transport proper- 
ties at low-frequency near the instability driven by com- 
pression. Including compression to the present effective 
medium approach will assess this view, and may provide 
a unified description of transport applicable to a broad 
range of amorphous solids. 

It is a pleasure to thank V. Vitelli for discussions stim- 
ulating this work, A. Liu, S. Nagel and N. Xu for discus- 
sions, and D. Huse and M. Muller for comments on the 
manuscript. 
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